本部分的重点包括:
- 线性回归的假定
- 探索性数据分析图和汇总统计数据
- 使用残差诊断回归假设
- 解释多元线性回归的参数和相关检验和区间
- 自助法置信区间的原理
\[\begin{aligned} y_i & = X_i \beta+\epsilon_i \\ & = \beta_1 X_{i1}+\cdots + \beta_k X_{ik} + \epsilon_i \end{aligned},\qquad for \quad i=1,\cdots , n\]
其中, \(\epsilon_i\) 相互独立,服从均值为0,标准差为 \(\sigma\) 正态分布
或者
\[y_i \sim N(X_i \beta, \sigma^2),\qquad for \quad i=1,\cdots , n\]
其中, \(X\) 是 \(n\times k\) 的预测变量矩阵,第 \(i\) 行为 \(X_{i}\), \(\beta\) 是长度为 \(k\) 的列向量
| 假定 | 表达式 | 假设不成立 | 解决方案 |
|---|---|---|---|
| 1 无完全共线性 | \(rank(X)=k,k<n\) | 无法识别回归系数 | 减少变量或增加样本 |
| 2 X外生 | \(E(X\epsilon )=0\) | 有偏,且不随N增加而改善 | 改变研究设计,因果推断技术 |
| 3 误差项零均值 | \(E(\epsilon)=0\) | 有偏,且不随N增加而改善 | 改变模型设定 |
| 4 线性关系 | \(y=\beta_{1}x_{1}+\beta_{2} x_{2}+\cdots+\epsilon\) | 模型设定错误 | 改变模型设定 |
| 5 误差项不相关(观测独立) | \(E(\epsilon_{i} \epsilon_{j})=0, i\neq j\) | 无偏但无效,标准误不正确 | 多层线性模型 |
| 6 同方差 | \(E(\epsilon^{'}\epsilon)=\sigma^{2}\boldsymbol{I}\) | 无偏但无效,标准误不正确 | 广义线性模型或稳健性标准误 |
| 7 误差项服从正态分布 | \(\epsilon\sim N(0,\sigma^{2})\) | 标准误不正确,但随N增加而改善 | 广义线性模型 |
Figure 1: 线性回归最小二乘法的假定
哪些研究设计可能导致违背回归假定?
如果违背正态性假定,那么需要广义线性模型;如果违背独立性假定,那么需要采用分层线性模型。
示例:母亲教育与子女认知能力数据集
library("foreign")
library("dplyr")
library("kableExtra")
library(gridExtra)
kidiq <- read.dta("kidiq.dta")
kidiq <- kidiq %>% mutate(hsfactor = ifelse(mom_hs==1, "高中及以上", "高中以下"))
table1 <- kidiq %>%
filter(row_number() < 6 | row_number() > 428)
kable(table1, booktabs=T,caption="示例数据概览") %>%
kable_styling(latex_options = "scale_down")
## Warning in attr(x, "align"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")
## Warning in attr(x, "align"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")
| kid_score | mom_hs | mom_iq | mom_work | mom_age | hsfactor | |
|---|---|---|---|---|---|---|
| 1 | 65 | 1 | 121.11753 | 4 | 27 | 高中及以上 | |
| 2 | 98 | 1 | 89.36188 | 4 | 25 | 高中及以上 | |
| 3 | 85 | 1 | 115.44316 | 4 | 27 | 高中及以上 | |
| 4 | 83 | 1 | 99.44964 | 3 | 25 | 高中及以上 | |
| 5 | 115 | 1 | 92.74571 | 4 | 27 | 高中及以上 | |
| 429 | 93 | 0 | 74.86073 | 2 | 25 | 高中以下 | |
| 430 | 94 | 0 | 84.87741 | 4 | 21 | 高中以下 | |
| 431 | 76 | 1 | 92.99039 | 4 | 23 | 高中及以上 | |
| 432 | 50 | 0 | 94.85971 | 2 | 24 | 高中以下 | |
| 433 | 88 | 1 | 96.85662 | 2 | 21 | 高中及以上 | |
| 434 | 70 | 1 | 91.25334 | 2 | 25 | 高中及以上 | |
有哪些变量类型?
Figure 2: 示例中尺度变量的直方图
这些变量的分布有什么特征?
library(GGally)
gg <- ggpairs(data = kidiq,
columns = c("hsfactor", "kid_score", "mom_iq", "mom_age"))
gg[4,1] <- gg[4,1] + geom_histogram( binwidth = 5)
gg[2,1] <- gg[2,1] + geom_histogram( binwidth = 10)
gg[3,1] <- gg[3,1] + geom_histogram( binwidth = 5)
gg
Figure 3: 探索变量间的关系
对角线上为单变量的分布,箱线图和直方图为按是否上过高中分类的其他尺度变量的分布,散点图为两个尺度变量的相关性,对角线上方为相关系数。
变量之间的关系是怎样的?
# Coded scatterplot
ggplot(kidiq, aes(x = mom_iq, y = kid_score, colour = hsfactor)) +
geom_point(aes(shape = hsfactor)) +
geom_smooth(aes(linetype = hsfactor), method = lm, se = FALSE)
Figure 4: 子女测试分数的线性趋势
通过探索性分析,可以尝试问这样的问题:
fit1 <- lm (kid_score ~ mom_iq, data = kidiq)
summary(fit1)
##
## Call:
## lm(formula = kid_score ~ mom_iq, data = kidiq)
##
## Residuals:
## Min 1Q Median 3Q Max
## -56.753 -12.074 2.217 11.710 47.691
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.79978 5.91741 4.36 1.63e-05 ***
## mom_iq 0.60997 0.05852 10.42 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.27 on 432 degrees of freedom
## Multiple R-squared: 0.201, Adjusted R-squared: 0.1991
## F-statistic: 108.6 on 1 and 432 DF, p-value: < 2.2e-16
但是注意:有时候在保持其他变量不变的情况下是不可能变动某个变量的,例如回归中同时包括IQ和 \(IQ^{2}\),或者有交互项 mom_hs*mom_i q.
回归模型诊断的常用方法,plot函数会自动绘制4张诊断图
par(mfrow=c(2,2))
plot(fit1)
# Fitted models for Model 2 and Model 2Q
ggplot(kidiq, aes(x = mom_iq, y = kid_score)) +
geom_point() +
stat_smooth(method = "lm", formula = y ~ x,
se = FALSE, linetype = 1) +
stat_smooth(method = "lm", formula = y ~ x + I(x^2),
se = FALSE, linetype = 2)
Figure 5: 比较线性和2次曲线拟合
kidiq <- kidiq %>% mutate(mom_iq2 = mom_iq*mom_iq)
fit2 <- lm (kid_score ~ mom_iq + mom_iq2, data = kidiq)
summary(fit2)
##
## Call:
## lm(formula = kid_score ~ mom_iq + mom_iq2, data = kidiq)
##
## Residuals:
## Min 1Q Median 3Q Max
## -54.824 -11.640 2.883 11.372 50.813
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -99.033675 37.301385 -2.655 0.008226 **
## mom_iq 3.076800 0.730291 4.213 3.07e-05 ***
## mom_iq2 -0.011917 0.003517 -3.389 0.000767 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.05 on 431 degrees of freedom
## Multiple R-squared: 0.2217, Adjusted R-squared: 0.2181
## F-statistic: 61.38 on 2 and 431 DF, p-value: < 2.2e-16
par(mfrow=c(2,2))
plot(fit2)
fit3 <- lm(kid_score ~ mom_hs, data = kidiq)
summary(fit3)
##
## Call:
## lm(formula = kid_score ~ mom_hs, data = kidiq)
##
## Residuals:
## Min 1Q Median 3Q Max
## -57.55 -13.32 2.68 14.68 58.45
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 77.548 2.059 37.670 < 2e-16 ***
## mom_hs 11.771 2.322 5.069 5.96e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 19.85 on 432 degrees of freedom
## Multiple R-squared: 0.05613, Adjusted R-squared: 0.05394
## F-statistic: 25.69 on 1 and 432 DF, p-value: 5.957e-07
此种情况下,回归系统和截距很容易解释。这个回归和两个独立样本t检验一致吗?需要满足等方差的假设吗?
fit4 <- lm(kid_score ~ mom_hs+ mom_iq, data = kidiq)
summary(fit4)
##
## Call:
## lm(formula = kid_score ~ mom_hs + mom_iq, data = kidiq)
##
## Residuals:
## Min 1Q Median 3Q Max
## -52.873 -12.663 2.404 11.356 49.545
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.73154 5.87521 4.380 1.49e-05 ***
## mom_hs 5.95012 2.21181 2.690 0.00742 **
## mom_iq 0.56391 0.06057 9.309 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.14 on 431 degrees of freedom
## Multiple R-squared: 0.2141, Adjusted R-squared: 0.2105
## F-statistic: 58.72 on 2 and 431 DF, p-value: < 2.2e-16
除了了解效应的显著性,还需要采用置信区间量化效应量的不确定性和量化模型预测值的不确定性。
confint(fit4)
## 2.5 % 97.5 %
## (Intercept) 14.1839148 37.2791615
## mom_hs 1.6028370 10.2973969
## mom_iq 0.4448487 0.6829634
new.data <- data.frame(mom_iq=190, mom_hs = 0)
predict(fit4, new = new.data, interval = "prediction")
## fit lwr upr
## 1 132.8737 95.18156 170.5658
如果回归模型违背了假设,自助法是一种更稳健的统计推断方法。其原理是原始数据样本能够代表总体,所以可以通过从原始样本数据中再抽样,获得多个子样本,从子样本的统计特征了解估计参数的不确定性。
自助法的步骤如下:
library(rsample)
library(tidyverse)
set.seed(666)
bootreg <- kidiq %>%
bootstraps(1000) %>%
pull(splits) %>%
map(\(df) lm(kid_score ~ mom_iq + mom_hs, data = df) ) %>%
map(\(mod) as.data.frame(t(as.matrix(coef(mod))))) %>%
list_rbind() %>% pivot_longer(everything(), names_to = "variable", values_to = "value")
bootregresult <- bootreg %>%
group_by(variable) %>%
summarise(across( where(is.numeric),
list(low = ~quantile(., probs = 0.025, na.rm = TRUE),
high = ~quantile(., probs = 0.975, na.rm = TRUE))))
bootregresult
# A tibble: 3 × 3
variable value_low value_high
<chr> <dbl> <dbl>
1 (Intercept) 14.1 37.5
2 mom_hs 1.51 10.5
3 mom_iq 0.439 0.681
ggplot(bootreg, aes(x = value)) +
geom_histogram(fill = "steelblue", color = "black", bins = 30) +
facet_wrap(~variable, scales = "free")
Figure 6: 回归系数的自助法分布
结果与正态分布理论假定下的结果相似。此外,还有多种其他计算自助法置信区间的方法,但是百分位法是应用最广泛的。
fit5 <- lm (kid_score ~ mom_hs + mom_iq - 1, data = kidiq)
summary(fit5)
##
## Call:
## lm(formula = kid_score ~ mom_hs + mom_iq - 1, data = kidiq)
##
## Residuals:
## Min 1Q Median 3Q Max
## -55.348 -10.796 2.187 12.789 50.838
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## mom_hs 5.99194 2.25786 2.654 0.00825 **
## mom_iq 0.81524 0.01979 41.189 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.51 on 432 degrees of freedom
## Multiple R-squared: 0.9571, Adjusted R-squared: 0.9569
## F-statistic: 4817 on 2 and 432 DF, p-value: < 2.2e-16
fit6 <- lm (kid_score ~ mom_hs + mom_iq + mom_hs:mom_iq, data = kidiq)
summary(fit6)
##
## Call:
## lm(formula = kid_score ~ mom_hs + mom_iq + mom_hs:mom_iq, data = kidiq)
##
## Residuals:
## Min 1Q Median 3Q Max
## -52.092 -11.332 2.066 11.663 43.880
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -11.4820 13.7580 -0.835 0.404422
## mom_hs 51.2682 15.3376 3.343 0.000902 ***
## mom_iq 0.9689 0.1483 6.531 1.84e-10 ***
## mom_hs:mom_iq -0.4843 0.1622 -2.985 0.002994 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.97 on 430 degrees of freedom
## Multiple R-squared: 0.2301, Adjusted R-squared: 0.2247
## F-statistic: 42.84 on 3 and 430 DF, p-value: < 2.2e-16
回到图4,交互项体现了图中的什么关系?
fit6 <- lm (kid_score ~ mom_hs * mom_iq, data = kidiq)
summary(fit6)
##
## Call:
## lm(formula = kid_score ~ mom_hs * mom_iq, data = kidiq)
##
## Residuals:
## Min 1Q Median 3Q Max
## -52.092 -11.332 2.066 11.663 43.880
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -11.4820 13.7580 -0.835 0.404422
## mom_hs 51.2682 15.3376 3.343 0.000902 ***
## mom_iq 0.9689 0.1483 6.531 1.84e-10 ***
## mom_hs:mom_iq -0.4843 0.1622 -2.985 0.002994 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.97 on 430 degrees of freedom
## Multiple R-squared: 0.2301, Adjusted R-squared: 0.2247
## F-statistic: 42.84 on 3 and 430 DF, p-value: < 2.2e-16
最终多元线性回归模型的典型特征包括:
虽然在报告研究结果时,一般要求选择一个合理的最终模型,但是值得注意:
1 研究者在得出结论时通常会检查并考虑一系列相关的模型
2 不同的研究者有时会针对同一组数据选择不同的模型作为“最终模型”。“最终模型”的选择取决于许多因素,例如主要研究问题、建模目的、简约性和拟合模型质量之间的权衡、基本假设等。建模决策不应自动化或完全基于统计检验;学科领域知识应始终在建模过程中发挥作用。要能够捍卫所选择的任何最终模型,但不应该感到有压力要找到唯一的“正确模型”,尽管大多数好的模型都会得出类似的结论。
在比较不同的模型构建时,可以使用的检验方法:
一个可能的最终模型如下:
fit7 <- lm (kid_score ~ mom_hs + mom_iq + mom_iq2, data = kidiq)
summary(fit7)
##
## Call:
## lm(formula = kid_score ~ mom_hs + mom_iq + mom_iq2, data = kidiq)
##
## Residuals:
## Min 1Q Median 3Q Max
## -54.553 -10.940 2.502 11.095 48.710
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -88.541848 37.390233 -2.368 0.018324 *
## mom_hs 5.107323 2.207015 2.314 0.021131 *
## mom_iq 2.828771 0.734492 3.851 0.000135 ***
## mom_iq2 -0.010910 0.003526 -3.094 0.002104 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.96 on 430 degrees of freedom
## Multiple R-squared: 0.2313, Adjusted R-squared: 0.2259
## F-statistic: 43.12 on 3 and 430 DF, p-value: < 2.2e-16
anova(fit7, fit4)
## Analysis of Variance Table
##
## Model 1: kid_score ~ mom_hs + mom_iq + mom_iq2
## Model 2: kid_score ~ mom_hs + mom_iq
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 430 138670
## 2 431 141757 -1 -3087 9.5723 0.002104 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
预期寿命的影响因素
life <- read.csv(file = "./life.csv")
life <- na.omit(life)
fit <- lm(formula = lifeexp ~ gdpcap + school + civlib + wartime,data = life)
summary(fit)
##
## Call:
## lm(formula = lifeexp ~ gdpcap + school + civlib + wartime, data = life)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.2871 -3.0705 -0.0368 4.6629 12.1285
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 52.9671697 3.0379356 17.435 < 2e-16 ***
## gdpcap 0.0003216 0.0002551 1.260 0.211
## school 2.3866392 0.4102266 5.818 9.01e-08 ***
## civlib -0.7522863 0.4861266 -1.548 0.125
## wartime 1.0936902 5.5354286 0.198 0.844
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.43 on 90 degrees of freedom
## Multiple R-squared: 0.7509, Adjusted R-squared: 0.7399
## F-statistic: 67.84 on 4 and 90 DF, p-value: < 2.2e-16
残差的理想模式:
x <- runif(1000,1,100)
error <- rnorm(1000)
y <- 10+ 3*x +error
fit0 <- lm(y~x)
plot(x=x,y=fit0$residuals)
存在异方差时的模式:
x <- runif(1000,1,100)
error <- rnorm(1000,mean=0,sd=1*x)
y <- 10+ 3*x +error
fit0 <- lm(y~x)
plot(x=x,y=fit0$residuals)
library(lmtest); library(car)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:purrr':
##
## some
## The following object is masked from 'package:dplyr':
##
## recode
vc <- hccm(fit, type = "hc1")
se.corrected <- sqrt(diag(vc))
coeftest(fit, vcov=vc)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 52.96716974 3.29289359 16.0853 < 2.2e-16 ***
## gdpcap 0.00032156 0.00026753 1.2020 0.2325
## school 2.38663924 0.42745500 5.5834 2.473e-07 ***
## civlib -0.75228629 0.50324260 -1.4949 0.1384
## wartime 1.09369023 5.30329592 0.2062 0.8371
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(fit$fitted.values,fit$residuals, main="life expectancy against fitted Y")
有异方差问题吗?还有什么问题?
残差与拟合值相关,意味着什么?
分别检查残差与各个协变量的关系
plot(life$civlib,fit$residuals)
plot(life$wartime,fit$residuals)
plot(life$school,fit$residuals)
plot(life$gdpcap,fit$residuals)
fit <- lm(formula = lifeexp ~ I(log(gdpcap)) + I(log(school)) + civlib + wartime, data=life)
summary(fit)
##
## Call:
## lm(formula = lifeexp ~ I(log(gdpcap)) + I(log(school)) + civlib +
## wartime, data = life)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.5038 -1.9584 0.2405 2.0338 11.0816
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.1893 5.4439 2.423 0.0174 *
## I(log(gdpcap)) 5.3730 0.6867 7.825 9.35e-12 ***
## I(log(school)) 6.0431 0.9040 6.685 1.88e-09 ***
## civlib -0.1843 0.3084 -0.598 0.5517
## wartime -0.2174 3.6894 -0.059 0.9531
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.627 on 90 degrees of freedom
## Multiple R-squared: 0.8889, Adjusted R-squared: 0.8839
## F-statistic: 180 on 4 and 90 DF, p-value: < 2.2e-16
plot(fit$fitted.values,fit$residuals, main="life expectancy against fitted Y")
plot(life$civlib,fit$residuals)
plot(life$wartime,fit$residuals)
plot(life$school,fit$residuals)
plot(life$gdpcap,fit$residuals)
观察残差图,可以看出相关模式已经消失了,但是还有异方差的可能(喇叭形)
比较稳健性标准误
library(orgutils)
fitsum <- summary(fit)
vc <- hccm(fit, type = "hc1")
fitout <- data.frame(est=fitsum$coefficients[,1], se=fitsum$coefficients[,2], robust_se=sqrt(diag(vc)))
toOrg(round(fitout,digits=2))
## | row.names | est | se | robust_se |
## |----------------+-------+------+-----------|
## | (Intercept) | 13.19 | 5.44 | 5.23 |
## | I(log(gdpcap)) | 5.37 | 0.69 | 0.65 |
## | I(log(school)) | 6.04 | 0.9 | 0.81 |
## | civlib | -0.18 | 0.31 | 0.28 |
## | wartime | -0.22 | 3.69 | 3.43 |
| 协变量 | \(R^2\) | \(\hat{\sigma}\) |
|---|---|---|
| GDP, School, Civlib, Wartime | 0.75 | 5.43 |
| log(GDP), log(School), Civlib, Wartime | 0.88 | 3.63 |
相同的数据,相同的响应变量
每个值是什么含义?哪个模型好一些?
样本外拟合优度
交叉验证
library(boot)
##
## Attaching package: 'boot'
## The following object is masked from 'package:car':
##
## logit
# 注意要采用广义线性回归函数
fitglm <- glm(lifeexp~gdpcap+school+civlib+wartime, data=life)
# 5重交叉验证
cv.err <- cv.glm(life,fitglm,K=5)
# 报告交叉验证的预测误差:第1个是原始值,第2个是调整值(与留一法比较)
cv.err$delta
## [1] 31.6745 31.2306
| 协变量 | \(R^2\) | \(\hat{\sigma}\) | 5重交叉验证误差 |
|---|---|---|---|
| GDP, School, Civlib, Wartime | 0.75 | 5.43 | 32.3 |
| log(GDP), log(School), Civlib, Wartime | 0.88 | 3.63 | 13.8 |
CV预测误差一般比回归标准误要高,为什么?
模型设定正确与否远比正确预测更重要,为什么?失去32年的预期寿命 vs. 失去14年预期寿命哪个更严重?
进一步改进模型
fit <- lm(formula = lifeexp ~ I(log(gdpcap)) + I(log(school)) + I(log(civlib)) + wartime + I(wartime^2), data=life)
summary(fit)
##
## Call:
## lm(formula = lifeexp ~ I(log(gdpcap)) + I(log(school)) + I(log(civlib)) +
## wartime + I(wartime^2), data = life)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.7938 -1.7693 0.0006 2.0123 10.5086
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.0854 5.4620 1.480 0.14233
## I(log(gdpcap)) 5.8587 0.6872 8.525 3.60e-13 ***
## I(log(school)) 6.0283 0.8606 7.005 4.53e-10 ***
## I(log(civlib)) 0.1266 0.8955 0.141 0.88790
## wartime 36.3329 12.8663 2.824 0.00585 **
## I(wartime^2) -110.8191 36.8656 -3.006 0.00344 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.482 on 89 degrees of freedom
## Multiple R-squared: 0.8987, Adjusted R-squared: 0.8931
## F-statistic: 158 on 5 and 89 DF, p-value: < 2.2e-16
fitglm <- glm(lifeexp~ I(log(gdpcap)) + I(log(school)) + I(log(civlib)) + wartime + I(wartime^2), data=life)
cv.err <- cv.glm(life,fitglm,K=5)
cv.err$delta
## [1] 12.95672 12.77216
| 协变量 | \(R^2\) | \(\hat{\sigma}\) | 5重交叉验证误差 |
|---|---|---|---|
| GDP, School, Civlib, Wartime | 0.75 | 5.43 | 32.3 |
| log(GDP), log(School), Civlib, Wartime | 0.88 | 3.63 | 13.8 |
| log(GDP), log(School),log(Civlib), Wartime, \(Wartime^2\) | 0.90 | 3.48 | 12.6 |
本部分知识重点:
投三次硬币决定来上课还是在宿舍睡觉,每出现一次正面代表上课多一票,每出现一次反面代表睡觉多一票。结果是前两次为正面、最后一次为反面。你可能会怀疑硬币是不是有问题,正面出现的概率真是1/2吗?可以计算正面的概率不同取值时,该事件的概率。
| \(\mathbf{y}\) | \(\hat{\theta}\) | \(\theta^{1s} \times (1-\theta)^{0s}\) | \(f_B (\mathbf{y} \vert \hat{\theta})\) |
|---|---|---|---|
| \({1,1,0}\) | \(0.00\) | \(0.00^2\times (1-0.00)^1\) | \(0.000\) |
| \({1,1,0}\) | \(0.10\) | \(0.10^2\times (1-0.10)^1\) | \(0.009\) |
| \({1,1,0}\) | \(0.20\) | \(0.20^2\times (1-0.20)^1\) | \(0.032\) |
| \({1,1,0}\) | \(0.30\) | \(0.30^2\times (1-0.30)^1\) | \(0.063\) |
| \({1,1,0}\) | \(0.40\) | \(0.40^2\times (1-0.40)^1\) | \(0.096\) |
| \({1,1,0}\) | \(0.50\) | \(0.50^2\times (1-0.50)^1\) | \(0.125\) |
| \({1,1,0}\) | \(0.60\) | \(0.60^2\times (1-0.60)^1\) | \(0.144\) |
| \({1,1,0}\) | \(0.67\) | \(0.67^2\times (1-0.67)^1\) | \(0.148\) |
| \({1,1,0}\) | \(0.70\) | \(0.70^2\times (1-0.70)^1\) | \(0.147\) |
| \({1,1,0}\) | \(0.80\) | \(0.80^2\times (1-0.80)^1\) | \(0.128\) |
| \({1,1,0}\) | \(0.90\) | \(0.90^2\times (1-0.90)^1\) | \(0.081\) |
| \({1,1,0}\) | \(1.00\) | \(1.00^2\times (1-1.00)^1\) | \(0.000\) |
我们将以分布参数的函数来描述观察数据的联合分布概率称之为似然函数(likelihood function),表示为\(L(\mathbf{y};\theta)\)。 \[L=\theta^2(1-\theta)^1\] \[log L=2log\theta + 1log(1-\theta)\] \[\frac{\partial log L}{\partial \theta}=\frac{2}{\theta}-\frac{1}{(1-\theta)}=0\] \[\hat \theta=2/3\] 最大化似然函数的\(\theta\)值是最大似然估计(maximum likelihood estimate, MLE)
最大似然推断并不是将数据视为随机的、参数视为固定的,而是将数据视为固定的,寻找什么样的参数最有可能生成这个数据,即将数据的联合分布视作某个分布密度或概率分布函数的参数的函数。
2012年189个国家的二氧化碳排放量与人均GDP(购买力平价),如果构建简单模型:\(Y=\beta_0+\beta_1 X +\epsilon\),其中Y为二氧化碳排放量的对数,X为人均GDP的对数。
library(ProfileLikelihood)
wdi$lgdppc<-log(wdi$gdp.pc.ppp)
xx <- profilelike.lm(formula = log(co2.kt)~1, data=wdi, profile.theta="lgdppc",
lo.theta=0.84, hi.theta=1.15, length=500)
with(xx,
plot(theta,profile.lik,las=1,lty=1,lwd=3,
type="l",pch=19,xlab=substitute(beta[1]),
ylab="likelihood",yaxt="n",bty="l",main="Least Squares as MLE",
xlim=c(0.85,1.1))
)
abline(v=coef(out)[2],col="gray50",lwd=3)
abline(h=max(xx$profile.lik),col="gray50",lwd=4)
# 建立自变量(含常数项)与因变量矩阵
x <- cbind(1,as.matrix(log(wdi$gdp.pc.ppp))) # 增加1列含1的常数项
y <- as.matrix(log(wdi$co2.kt))
K <- ncol(x); n <- nrow(x) # 观测数n和变量数K
# 定义似然函数,可以选择多种参数化方式,此处采用logL的完整形式
loglik.my <- function(par,X,Y) {
Y <- as.vector(y)
X <- as.matrix(x)
xbeta <- X%*%par[1:K]
sigma <- sqrt(sum(((X[,2]-mean(X[,2]))^2)/(n-K))) # 假定标准误已知,有多种设定形式
sum(-(1/2)*log(2*pi)-(1/2)*log(sigma^2)-(1/(2*sigma^2))*(y-xbeta)^2) # 对数似然函数,加负号变为最小化
}
# 将似然函数传递给最优化函数,提供初始值,选择算法,设定迭代次数等
mle.fit <- optim(c(5,5),loglik.my, method = "BFGS", control =
list(trace=TRUE,maxit=10000,fnscale = -1),hessian = TRUE)
## initial value 150490.736339
## final value 375.927547
## converged
if(mle.fit$convergence!=0)
print("MDW WARNING: Convergence Problems; Try again!")
# 计算标准诊断量
stderrors<- sqrt(diag(-solve(mle.fit$hessian)))
z<-mle.fit$par/stderrors
p.z <- 2* (1 - pnorm(abs(z)))
out.table <- data.frame(Est=mle.fit$par, SE=stderrors, Z=z, pval=p.z)
round(out.table, 2)
## Est SE Z pval
## 1 2.69 0.67 4.02 0.00
## 2 0.18 0.07 2.48 0.01
线性回归模型假定同方差,即对于所有的观测都有 \(y_i \sim N(\mu_i , \sigma^2 )\),但是如果误差项是异方差的,即 \(y_i \sim N(\mu_i , \sigma_i^2 )\),则会导致\(se(\hat{\beta})\)是有偏的,\(\hat{\beta}\)的估计是无效的。
一般采用稳健性标准误来替代“有问题”的标准误,但是,异方差被看做是数据存在的一个“问题”,仅仅是因为我们假定它不应该在“正常”的数据中存在。真实的可能是我们的线性回归假定框定了我们解决问题的范围。类比一下,为什么我们不会说“异均值”是一个问题,是因为线性回归的框架中允许了因变量均值随自变量变动(这也正是我们建模的基础),那如果我们能把\(\sigma_i^2\) 也纳入到模型中,即对因变量方差建模,也就能更好的利用自变量解释解释因变量的均值和方差两方面的变化。
异方差暗含了自变量与因变量方差的关系,这可能正是我们想研究的内容。
异方差正态模型的最大似然估计
异方差最大似然估计的推导:
假想样本量为2000的异方差数据来自真实模型:
\[y_i \sim N (\mu_i , \sigma_i^2 )\]
\[\mu_i = 5 + 10x_i \]
\[\sigma_i^2 = exp(1+3x_i ) \]
分别采用线性回归与异方差的最大似然估计去拟合
set.seed(1234)
obs <- 2000
x <- runif(obs)
mu <- 5+10*x
sigma2 <- exp(1+3*x)
y <- rnorm(obs, mu, sqrt(sigma2))
# 线性回归拟合
lm.fit <- lm(y~x)
summary(lm.fit)
##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.0880 -2.2792 0.0612 2.2686 19.9045
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.0546 0.1839 27.49 <2e-16 ***
## x 10.0479 0.3206 31.34 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.108 on 1998 degrees of freedom
## Multiple R-squared: 0.3296, Adjusted R-squared: 0.3292
## F-statistic: 982.2 on 1 and 1998 DF, p-value: < 2.2e-16
# 最大似然法拟合
# 建立两个系统部分的自变量
xcovariates <- x
zcovariates <- x
# 设置参数初始值beta0,beta1,gamma0,gamma1
stval <- c(0,0,0,0)
# 设定似然函数
llk.hetnormlin <- function(param,y,x,z) {
x <- as.matrix(x)
z <- as.matrix(z)
x <- cbind(1,x)
z <- cbind(1,z)
b <- param[ 1 : ncol(x) ]
g <- param[ (ncol(x)+1) : (ncol(x) + ncol(z)) ]
xb <- x%*%b
s2 <- exp(z%*%g)
sum(0.5*(log(s2)+(y-xb)^2/s2)) # optim是最小化函数,所以公式要乘以-1
}
# 运行优化函数得到拟合结果
mle.fit <- optim(stval,llk.hetnormlin,method="BFGS",hessian=T,y=y,x=xcovariates,z=zcovariates)
# 提取估计值
pe <- mle.fit$par # 参数的点估计
vc <- solve(mle.fit$hessian) # 协方差矩阵
se <- sqrt(diag(vc)) # 标准误
z <- mle.fit$par/se # z值
p.z <- 2* (1 - pnorm(abs(z))) # p值
out.table <- data.frame(Coefs=c("beta0","beta1","gamma0","gamma1"), Est=mle.fit$par, SE=se, Z=z, pval=p.z)
print(out.table, digits=2)
## Coefs Est SE Z pval
## 1 beta0 5.05 0.102 50 0
## 2 beta1 10.05 0.278 36 0
## 3 gamma0 0.99 0.064 15 0
## 4 gamma1 3.01 0.112 27 0
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select